PTC Liverpool

1. Install the required packages

library(cyclestreets)
library(mapview)
library(pct)
library(sf)

## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3

library(stats19)

## Data provided under OGL v3.0. Cite the source and link to:
## www.nationalarchives.gov.uk/doc/open-government-licence/version/3/

library(stplanr)
library(tidyverse)

## -- Attaching packages ------------------------------------------------------------------------------- tidyverse 1.2.1 --

## v ggplot2 3.2.1     v purrr   0.3.2
## v tibble  2.1.3     v dplyr   0.8.3
## v tidyr   1.0.0     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.4.0

## -- Conflicts ---------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(tmap)
library(tmaptools)
library(dplyr)
library(leaflet)
library(classInt)
library(ggplot2)

2. Download the pct desire lines for Liverpool via the PCT

region_name = "liverpool-city-region"
max_distance = 7
lines_test = get_pct_lines(region_name)
zones_all = get_pct_zones(region_name)
L_original_lsoa = get_pct_lines("liverpool-city-region", geography = "lsoa")

Liverpool_filter = filter(L_original_lsoa, lad_name1 %in% c("Liverpool")) # filter out Liverpool only

3. Create bicycle and car dependent interative maps

tmap_mode("view")

## tmap mode set to interactive viewing

#3a Bicycle Dependent Desire Lines
bicycle = Liverpool_filter %>%
  mutate('Percentage Cycling' = (bicycle) / all * 100) %>%
  filter(e_dist_km < max_distance)
tm_shape(bicycle) +
  tm_lines("Percentage Cycling", palette = "RdYlBu", lwd = "bicycle", scale = 9)

## Legend for line widths not available in view mode.

#3b Car Dependent Desire Lines
car_dependent = Liverpool_filter %>% 
  mutate(`Percent Drive` = (car_driver) / all * 100) %>% 
  filter(e_dist_km < max_distance)
tm_shape(car_dependent) +
  tm_lines("Percent Drive", palette = "-RdYlBu", lwd = "all", scale = 9)+
  tm_layout(legend.text.size = 1.2)

## Legend for line widths not available in view mode.

#note: these desire lines can be created for any travel mode, adapt the code above to show desired mode

4. Extract the desire lines that fall within the top 10 cycling trips

Top10_Liverpool = Liverpool_filter %>% top_n(10, bicycle)

#filter the columns
Top10_Liverpool_data = Top10_Liverpool %>% select(id, geo_code1, geo_code2, geo_name1, geo_name2, lad11cd1, lad11cd2, lad_name1, lad_name2, all, bicycle, foot, car_driver, car_passenger, motorbike, train_tube, bus, taxi_other, govtarget_slc, dutch_slc)

5. Create routes using the Cycle Streets API

#you will require a cyclestreets API to run this code

#code = Top10_Liverpool_route = stplanr::line2route(Top10_Liverpool_data, route_fun = stplanr::route_cyclestreet)

Top10_Liverpool_route <- st_read("C:/Users/Holly.Mizser-Jones/OneDrive - Arup/UCL/CASA005/Assessment/Data/Top10_liverpool_route.shp")

## Reading layer `Top10_liverpool_route' from data source `C:\Users\Holly.Mizser-Jones\OneDrive - Arup\UCL\CASA005\Assessment\Data\Top10_liverpool_route.shp' using driver `ESRI Shapefile'
## Simple feature collection with 11 features and 20 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: -2.99466 ymin: 53.34071 xmax: -2.83573 ymax: 53.41017
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

#join the Government Target, Go Dutch and Existing Bicycle trips to the raw route data
Top10_Liverpool_route$govtarget = Top10_Liverpool_data$govtarget_slc
Top10_Liverpool_route$dutch = Top10_Liverpool_data$dutch_slc
Top10_Liverpool_route$bicycle <- Top10_Liverpool_data$bicycle

6. Create graph showing the difference by distance and hilliness per scenario

#map graph showing the variation in hilliness / distance between the Government Target and Go Dutch scenarios

#difference between gradient and distance

distances = 1:20
hilliness = 0:5
uptake_df = data.frame(
  distances = rep(distances, 6),
  hilliness = rep(hilliness, each = 20)
)
p_govtarget = uptake_pct_govtarget(
  distance = uptake_df$distances,
  gradient = uptake_df$hilliness
)
p_godutch = uptake_pct_godutch(
  distance = uptake_df$distances,
  gradient = uptake_df$hilliness
)
uptake_df = rbind(
  cbind(uptake_df, scenario = "govtarget", pcycle = p_govtarget),
  cbind(uptake_df, scenario = "godutch", pcycle = p_godutch)
)
library(ggplot2)
ggplot(uptake_df) +
  geom_line(aes(
    distances,
    pcycle,
    linetype = scenario,
    colour = as.character(hilliness)
  )) +
  scale_color_discrete("Gradient (%)")+
  labs(x = "Distance (km)", y = "Propensity to Cycle (Multiplier)")

7. Health Data

#read in LSOA shapefile for UK with IMD 2019 data joined
Liverpool_lsoa_shp <- st_read("C:/Users/Holly.Mizser-Jones/OneDrive - Arup/UCL/CASA005/Assessment/Shp/Indices_of_Multiple_Deprivation_IMD_2019/Indices_of_Multiple_Deprivation_IMD_2019.shp")

## Reading layer `Indices_of_Multiple_Deprivation_IMD_2019' from data source `C:\Users\Holly.Mizser-Jones\OneDrive - Arup\UCL\CASA005\Assessment\Shp\Indices_of_Multiple_Deprivation_IMD_2019\Indices_of_Multiple_Deprivation_IMD_2019.shp' using driver `ESRI Shapefile'
## Simple feature collection with 32844 features and 66 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -6.418524 ymin: 49.86474 xmax: 1.762942 ymax: 55.81107
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

#filter for Liverpool
Liverpool_filter_shp = filter(Liverpool_lsoa_shp, LADnm %in% c("Liverpool"))

#visualise the routes on top of the IMD polygon data
library(leaflet)

pal <- colorFactor(get_brewer_pal("-YlGnBu", n = 10), domain = Liverpool_filter_shp$IMD_Decile, n = 10)

pal2 <- colorFactor(get_brewer_pal("-YlGnBu", n = 10), domain = Liverpool_filter_shp$HDDDec, n = 10)

IMDandroutes<- leaflet() %>%
  # add basemap options
  addProviderTiles(providers$Esri.WorldGrayCanvas, group = "Esri Grey") %>%
  addTiles(group = "OSM") %>%
  
  #add our Borough polygons, linking to the tables we just made
  addPolygons(data=Liverpool_filter_shp,
              color="white",
              weight = 2,
              opacity = 1,
              dashArray = "3",
              fillOpacity = 0.6,
              fillColor = ~pal(Liverpool_filter_shp$IMD_Decile),
              group = "Index of Multiple Deprivation")%>%
  
  # add a legend for boroughs
  addLegend(pal = pal, 
            values = Liverpool_filter_shp$IMD_Decile,
            group=c("Index of Multiple Deprivation"), 
            title ="IMD Decile",
            position ="bottomleft")%>%
  
  #add our ward polygons, linking to the tables we just made
  addPolygons(data=Liverpool_filter_shp,
              color="white",
              weight = 2,
              opacity = 1,
              dashArray = "3",
              fillOpacity = 0.6,
              fillColor = ~pal2(Liverpool_filter_shp$HDDDec),
              group = "Health Deprivation and Disability Decile")%>%
  
  # add a legend for wards
  addLegend(pal = pal2, 
            values = Liverpool_filter_shp$HDDDec,
            group=c("Health Deprivation and Disability Decile"), 
            position ="bottomleft",
            title ="HDD Decile")%>%
 
   #add routes
  addPolylines(data=Top10_Liverpool_route, color = "black", 
               weight = 2.5, opacity = 1, group = c("Top 10 Cycle Routes in Liverpool")) %>%

  # specify layers control
  addLayersControl(
    baseGroups = c("Esri Grey Canvas", "OSM"),
    overlayGroups = c("Index of Multiple Deprivation", "Health Deprivation and Disability Decile","Top 10 Cycle Routes in Liverpool"),
    options = layersControlOptions(collapsed = FALSE))%>%
  hideGroup(c("Index of Multiple Deprivation", "Health, Deprivation and Disability Decile", "Top 10 Cycle Routes in Liverpool"))

# show us the map
IMDandroutes

7a Health Intersect

#intersect the Top Routes with the IMD polygon data
IMD_route_intersect <- st_join(Top10_Liverpool_route, Liverpool_filter_shp, join = st_intersects)

## although coordinates are longitude/latitude, st_intersects assumes that they are planar

#export this to csv and use Microsoft Excel to calculate the average for each route
write.csv(IMD_route_intersect,"C:/Users/Holly.Mizser-Jones/Documents/UCL/CASA005/Assessment//IMD_route_average.csv", row.names = FALSE)

#read the IMD average back into R
IMD_route_average <- read.csv("C:/Users/Holly.Mizser-Jones/Documents/UCL/CASA005/Assessment/IMD_average.csv")

#merge with the route data 
route_merge <-merge (Top10_Liverpool_route, 
                     IMD_route_average, 
                     by.x="length", 
                     by.y="length",
                     no.dups = TRUE)

7b Calculate the change in Marginal METS (Goodman et al., 2019)

#the change in mMETS was calculated in excel, the below code joins the resultant increases to the routes

#link MET data to Top_10_routes

mMET <- read.csv("C:/Users/Holly.Mizser-Jones/OneDrive - Arup/UCL/CASA005/Assessment/liverpool_mMET.csv")

#merge mMET with Routes
route_mMET <-merge (route_merge, 
                     mMET, 
                     by.x="Route", 
                     by.y="Route",
                     no.dups = TRUE)